{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chris Holden (ceholden@gmail.com) - https://github.com/ceholden\n",
    "\n",
    "Chapter 2: Your first remote sensing vegetation index\n",
    "====================================================="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "Now that we can read our data into the computer, let's calculate some vegetation indices.\n",
    "\n",
    "The [Normalized Difference Vegetation Index (NDVI)](http://en.wikipedia.org/wiki/NDVI) is so ubiquitous that it even has a Wikipedia entry. If you're here for learning how to do remote sensing image processing using GDAL and Python, I suspect you don't need any introduction to this section. If you need a refresher, please visit the Wikipedia URL for [NDVI](http://en.wikipedia.org/wiki/NDVI)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This chapter will be very straightfoward. We've already seen how we can read our imagery into a NumPy array -- this chapter will simply extend these tools by showing how to do simple calculations on NumPy objects.\n",
    "\n",
    "Let's bring up our previous code for opening our image and reading in the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Red band mean: 589.379808\n",
      "NIR band mean: 3442.297712\n"
     ]
    }
   ],
   "source": [
    "# Import the Python 3 print function\n",
    "from __future__ import print_function\n",
    "\n",
    "# Import the \"gdal\" and \"gdal_array\" submodules from within the \"osgeo\" module\n",
    "from osgeo import gdal\n",
    "from osgeo import gdal_array\n",
    "\n",
    "# Import the NumPy module\n",
    "import numpy as np\n",
    "\n",
    "# Open a GDAL dataset\n",
    "dataset = gdal.Open('../../example/LE70220491999322EDC01_stack.gtif', gdal.GA_ReadOnly)\n",
    "\n",
    "# Allocate our array using the first band's datatype\n",
    "image_datatype = dataset.GetRasterBand(1).DataType\n",
    "\n",
    "image = np.zeros((dataset.RasterYSize, dataset.RasterXSize, dataset.RasterCount),\n",
    "                 dtype=gdal_array.GDALTypeCodeToNumericTypeCode(image_datatype))\n",
    "\n",
    "# Loop over all bands in dataset\n",
    "for b in range(dataset.RasterCount):\n",
    "    # Remember, GDAL index is on 1, but Python is on 0 -- so we add 1 for our GDAL calls\n",
    "    band = dataset.GetRasterBand(b + 1)\n",
    "    \n",
    "    # Read in the band's data into the third dimension of our array\n",
    "    image[:, :, b] = band.ReadAsArray()\n",
    "    \n",
    "\n",
    "print('Red band mean: {r}'.format(r=image[:, :, 2].mean()))\n",
    "print('NIR band mean: {nir}'.format(nir=image[:, :, 3].mean()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even from simple mean statistics over the entire image we can see the contrast between the red and the near-infared (NIR) bands.\n",
    "\n",
    "## NDVI\n",
    "\n",
    "To calculate NDVI, we can simply use standard arithmetic operators in Python because these operations in NumPy are vectorized. Just like MATLAB, R, and other higher level languages, **NEVER** loop over a NumPy array unless you can avoid it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.71390828  0.71079741  0.69352291 ...,  0.79392185  0.81408451\n",
      "   0.79165379]\n",
      " [ 0.68064263  0.6787194   0.6643924  ...,  0.81387182  0.79880597\n",
      "   0.77389811]\n",
      " [ 0.66904762  0.67268446  0.66332892 ...,  0.78495923  0.78278801\n",
      "   0.81253291]\n",
      " ..., \n",
      " [ 0.68301262  0.68593651  0.67145614 ...,  0.81065089  0.78050922\n",
      "   0.76519266]\n",
      " [ 0.67341718  0.6622986   0.65331611 ...,  0.80436681  0.77483099  0.75      ]\n",
      " [ 0.63973799  0.62396514  0.66731813 ...,  0.7094648   0.70005244\n",
      "   0.74574523]]\n",
      "0.904601300891\n"
     ]
    }
   ],
   "source": [
    "b_red = 2\n",
    "b_nir = 3\n",
    "\n",
    "ndvi = (image[:, :, b_nir] - image[:, :, b_red]) / (image[:, :, b_red] + image[:, :, b_nir])\n",
    "\n",
    "print(ndvi)\n",
    "print(ndvi.max())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Note: Python 2\n",
    "\n",
    "In Python 2 an integer divided by an integer produces an integer, even if the division would have produced a float point number. Python 3 changed this behavior, but if we run the NDVI calculation with Python 2 we would end up with all of our NDVI values equal to 0 because our input image is an integer datatype (int16). See [documentation for division in NumPy for more information](http://docs.scipy.org/doc/numpy/reference/generated/numpy.divide.html).\n",
    "\n",
    "While we don't necessarily need to change anything in Python 3, it is generally useful to be explicit with the datatypes involved in our calculations for the sake of clarity. Additionally, we generally also want code written using Python 3 to work with Python 2.\n",
    "\n",
    "In order to ensure we perform our calculation using floating point numbers, we can cast the demoninator or numerator of the calculation to a float:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NDVI matrix: \n",
      "[[ 0.71390828  0.71079741  0.69352291 ...,  0.79392185  0.81408451\n",
      "   0.79165379]\n",
      " [ 0.68064263  0.6787194   0.6643924  ...,  0.81387182  0.79880597\n",
      "   0.77389811]\n",
      " [ 0.66904762  0.67268446  0.66332892 ...,  0.78495923  0.78278801\n",
      "   0.81253291]\n",
      " ..., \n",
      " [ 0.68301262  0.68593651  0.67145614 ...,  0.81065089  0.78050922\n",
      "   0.76519266]\n",
      " [ 0.67341718  0.6622986   0.65331611 ...,  0.80436681  0.77483099  0.75      ]\n",
      " [ 0.63973799  0.62396514  0.66731813 ...,  0.7094648   0.70005244\n",
      "   0.74574523]]\n",
      "\n",
      "Max NDVI: 0.9046013008913515\n",
      "Mean NDVI: 0.7088133953809207\n",
      "Median NDVI: 0.7319195214790647\n",
      "Min NDVI: 0.09470304975922954\n"
     ]
    }
   ],
   "source": [
    "ndvi = (image[:, :, b_nir] - image[:, :, b_red]) / \\\n",
    "        (image[:, :, b_nir] + image[:, :, b_red]).astype(np.float64)\n",
    "\n",
    "print('NDVI matrix: ')\n",
    "print(ndvi)\n",
    "\n",
    "print('\\nMax NDVI: {m}'.format(m=ndvi.max()))\n",
    "print('Mean NDVI: {m}'.format(m=ndvi.mean()))\n",
    "print('Median NDVI: {m}'.format(m=np.median(ndvi)))\n",
    "print('Min NDVI: {m}'.format(m=ndvi.min()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks correct.\n",
    "\n",
    "Speaking of *looking correct*, the next chapter (link to [webpage](chapter_3_visualization.html) or [Notebook](chapter_3_visualization.ipynb)) will demonstrate how you can visualize your results using actual plots!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
